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^ . Progress in numerical relativity has been hindered for 30 years because of the 

^^ ! difficulties of avoiding spacetinie singularities in numerical evolution. We propose a 

Q^ I scheme which excises a region inside an apparent horizon containing the singularity. 

Two major ingredients of the scheme are the use of a horizon- locking coordinate and a 
finite differencing which respects the causal structure of the spacetime. Encouraging 
vQ . results of the scheme in the spherical collapse case are given. 
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I. INTRODUCTION 

Numerical studies of the Einstein equations promise to deepen dramatically our knowl- 
edge of general relativity, astrophysics, and cosmology. Shapiro and Teukolsky |Q suggested 
that the "Holy Grail of numerical relativity" is a code that (i) avoids singularities, (ii) han- 
dles black holes, (iii) maintains high accuracy, and (iv) is capable of running forever. These 
four difficulties in numerical relativity are closely tied to each other. Black holes, and within 
them spacetime singularities, may be formed after some evolution, even if none exists ini- 
tially. The presence of such objects implies that the dynamical range of the calculation is 
large, with very different length and time scales involved. This makes it difficult to maintain 
accuracy or even to keep the code from crashing in numerical evolution. 

The traditional way to deal with this problem is to use the coordinate degrees of freedom 
to avoid the "extreme" regions. With "the many fingers of time" in relativity, one may 
evolve other regions in space without evolving the region in which a singularity is about to 
form. Many different types of singularity-avoiding slicings have been proposed and studied 
in detail (see, e.g. [§-§]). This idea of using the freedom in slicing the spacetime to avoid 
singularities is ingenious but not perfect. In the vicinity of the singularity these slicings 
inevitably contain a region of abrupt change near the horizon, and a region in which the 
constant time slices dip back deep into the past in some sense. Depending on the details 
of the choices of the spacetime coordinates, the code will sooner or later crash due to these 
pathological properties of the slicing. The problem can appear as the development of spikes 
in the spatial metric functions , the steepening of spatial gradients |l[ , "grid stretching" |l| 
or large coordinate shift P] on the black hole throat, etc. 

Hence it is important to investigate other ways to handle singularities and black holes in 
numerical relativity. Cosmic censorship suggests that in physical situations, singularities are 
hidden inside black hole horizons. With the problematic region of numerical evolution mostly 
inside the horizon, it is tempting to cut away this region from the numerical calculation by 
imposing a boundary condition on or slightly inside the horizon. To an outside observer no 
information is lost since the region cut away is unobservable. There is then no singularity 
to crash the code, the observable region can be evolved, the dynamic range is drastically 
reduced so accuracy is easier to maintain, and there is in principle no physical reason that 
the code cannot run forever. 

Using a horizon boundary condition by itself is not a new idea m^-^. However, it is non- 
trivial to implement a horizon boundary condition in dynamical evolution [^j . The boundary 
condition to be imposed on a black hole horizon, which is a one way membrane ||T0[ , should 



presumably be some sort of out-going (into the hole) boundary condition. However, except 
for the case of linear, non-dispersive fields propagating in a fiat spacetime, we are not aware 



of any satisfactory out-going wave boundary condition in numerical calculation ||Tl[, let 
alone for waves in relativity, which can be non-linear, dispersive, and have tails and other 
complications [|12] . 



In this Letter we demonstrate that the idea of a horizon boundary condition in numerical 
relativity can indeed be realized. In the next section, we discuss the two major ingredients 
for our successful implementation of it: (i) a "horizon locking coordinate" (HLC) system 
which ties the spatial coordinates to the spatial geometry, (ii) a finite differencing scheme 
that we call causal differencing (CD) which respects the causal structure of the spacetime. 



CD is not only essential for the stability of the code using the HLC, but also eliminates the 
need of explicitly imposing boundary conditions on the horizon. It is, in a sense, the horizon 
boundary condition without a boundary condition: since the horizon is a one way membrane, 
the quantities on the horizon can be affected only by quantities outside but not inside the 
horizon. Hence, in CD, all quantities on the boundary can be updated explicitly in terms 
of known quantities residing on or outside the boundary, provided the boundary is inside 
the horizon. There is no need to impose boundary conditions to account for information 
not covered by the numerical evolution. Such a horizon boundary condition is particularly 
desirable since it is general to all kinds of source terms in the Einstein equations. 

A. Horizon Locking Coordinate 

We use the numerical construction of the Schwarzschild spacetime to illustrate the idea 
of locking the horizon. As shown in 0, the numerical construction of a Schwarzschild 
spacetime is non-trivial. In all nine differencing schemes and for all choices of coordinate 
conditions studied in |^, problems arise after an evolution in time of about lOOM with 
reasonable choices of grid parameters. To facilitate comparison with the results of 0, we 
write the line element in the same form: 

ds"^ = {-a^ + ilj^Af3^)dt^ + 2^^Al3dtd'q + ^^[Adrf + BTf{de^ + sin^O rf^^)] (1) 

Such a line element is easily generalized to one which is suitable for numerical study of 
axisymmetric spacetimes |T3|, and it includes both the radial gauge |T^ and the quasi- 
isotropic or isothermal gauge [p], |l5|] . 



In p[, the initial data used in evolving the Schwarzschild geometry are determined by 



time symmetry and conformal flatness [|l^, that is, the initial slice is an Einstein- Rosen 
bridge. They find that the most successful code is one which uses the MacCormack or 
Brailovskaya differencing scheme, maximal slicing [0 and zero shift. This code (denoted by 
BHS) is considered to be very accurate, but as for all codes designed so far to handle black 
holes, it will develop difficulties at late times. In Fig.l, the dashed line labeled "Free Horizon" 
gives the coordinate position of the apparent horizon (AH) vs. time. We see that the AH 
is rapidly growing in radius. The solid lines give the value of the radial metric component 
A every t = lOM. A spike is rapidly developing near the AH, which eventually causes the 
code to become inaccurate and crash. The inaccuracy generated by the sharp spike shows 
up both in the violation of the Hamiltonian constraint, (dashed line with asterisks in Fig. 
2) and the drifting of the mass of the black hole from its initial value 2 |T7 . 



A natural way to avoid this development of a sharp peak is to tie the coordinate grid 
points to some geometric structure. To achieve this, (i)We make the coordinate position of 
the AH constant in time after it grows to a certain finite radius. This determines the shift (3 
at the AH. (ii)After the horizon is locked, all grid points are tied to it by requiring the radial 
metric function A{t, rf) to be a constant in time. This determines the shift elsewhere, (iii) 
We drop most grid points from the calculation inside the AH, leaving only a small buffer 
zone. 



The results using the HLC are shown in the Figs. [fq]. Fig. 1 shows that after the shift 



is phased in, the AH is locked and instead of developing a spike as in BHS or other codes. 



the radial metric function A is 0(1) everywhere and absolutely unchanged over time. In 
the HLC with the singular region cut away, there is no need to use any singularity- avoiding 
time slicing. To demonstrate this freedom, we used a lapse that is constant in time after 



t = 5M (but a > 0.3 everywhere), when the AH is securely locked |T^. Without the sharp 
spike in the metric function. Fig. 2 shows that outside the AH the Hamiltonian constraint 
is satisfied to more than an order of magnitude better than in BHS by t = lOOM, and that 
the mass of the black hole remains essentially constant for all time as it should be in this 
vacuum case. 

B. Causal Differencing 

One consequence of introducing a shift vector in the HLC is that inside the AH the future 
light cone is tilted inward towards smaller rj. This feature is essential for our implementation 
of the AH condition. It allows us to do without a specific boundary condition for each type 
of ingoing wave since inside the AH, data at a particular grid point depend explicitly only 
on past data from grid points at equal or larger coordinate points. The existence of a shift 
calls for a finite differencing scheme different from the usual one. Here we illustrate this 
with a simple example. 

Consider a scalar field in a 1+1 flat space: ds"^ = —df'^ + dx''^ and df,ip = d'^,(f). Suppose 
we break it up for numerical evolution as [^ dt'(f) = d^'ir and dt'ir = dx'4>. If we introduce 
a constant shift f3 by the coordinate change t = t' and x = x' — (3t', then the wave equation 
becomes dt4> = dxTi + Pdx4> and dtn = dx(p + PdxTi- The question is how to write down the 
finite difference version of these two equations with a shift. One might naively use, say, the 
usual leapfrog scheme: 

2At 2Ax '' 2Ax ^ ^ 

and likewise for the dfTT equation. However, a von Neumann stability analysis shows that for 
any given C = At/ Ax, the scheme is unstable for a large enough j3. For example, numerical 
experiments show that for C = 1/3, the instability quickly takes over with jSC = 0.67. 

To obtain a stable finite differencing, we do the following: (i) Return to the unshifted wave 
equations in the primed coordinate, which have untilted light cones. The finite differencing 
can be written as usual, (ii) Transform directly the finite difference equation into the 
unprimed coordinates. In leapfrog scheme, this leads to 

2At 2 Ax ^^ 2Ax' ^ ' 

where Ax' = (3 At. Notice that on the RHS, the spatial derivative is centered at j + (3C, 
which is the center (on the n*^ slice) of the causal dependence of the point (n+l,j) where is 
to be updated. This guarantees that the region of causal dependence is completely covered 
whenever C < 1. 

The above procedure of obtaining a finite difference scheme which observes the causal 
structure in the presence of a shift vector can be applied to any difference scheme, just 
as for the leapfrog example here. The underlying idea of CD is similar to the "up- wind" 



differencing scheme used in hydrodynamic calculations ||2T| p2[. Several variations on the 



form of the CD operator will be discussed in a future paper. 

II. THE APPARENT HORIZON SCHEME IN SPHERICAL COLLAPSE 

In this section we demonstrate these ideas with a spherical collapse code. A self- 
gravitating massless Klein-Gordon scalar field in an initial time-symmetric gaussian distri- 
bution is added to the black hole. Theoretical and computational details will be presented 
in a future paper. Here we summarize results for a typical case. 

In Fig. 2b solid (dashed) lines represent the AH mass M^h as a function of time for 
the cases with (without) the AH scheme. The ADM mass of the entire spacetime for the 
case presented here was Madm = 2.62, while initially Mah = 2.08. With the gaussian 
distribution of the matter initially centered at rj = 2.5 and a width of Srj = 0.5, the hole 
has absorbed nearly all ingoing radiation by t = 15. We see that with the AH scheme, the 
solid line levels off after that time as expected, with a mass of 2.33. We see no reflection of 
radiation from the horizon. The dashed line shows results with BHS where Mah continues 
to drift upward. In fact, because of the large spike which develops in the metric function 



A, the system becomes unstable at about t = 80M . The comparison of the accuracy ||2^ 
between the AH and BHS schemes is also shown in Fig. 2a. Again we see that the violation 
of the hamiltonian constraint is increasing in the BHS case, while it is essentially constant 
over time with the AH scheme, suggesting that the code is capable of running for a long 
time. One final note is that the AH scheme is also potentially much faster than one using 
maximal slicing, which involves solving an elliptic equation. 

III. CONCLUSIONS 

Progress in numerical relativity has been hindered for 30 years because of the difficulties 
of avoiding spacetime singularities in the calculations. We have presented working exam- 
ples of how an apparent horizon boundary scheme can help circumvent these difficulties in 
dynamic spacetimes. 

There are a number of issues which must be addressed in future work. For example, 
in some codes the constraint equations are solved explicitly for some components of the 
extrinsic curvature or the metric during the evolution. These elliptic equations would require 
boundary conditions which may be difficult to formulate in an AH scheme. However, in a 
free evolution scheme such issues do not arise. Also, it is possible for a time slice to hit a 
singularity before an AH is formed (regardless of questions about Cosmic Censorship) [Q. 
This potential problem can be handled by using a singularity avoiding slicing condition until 
the horizon is formed and locked in place, as we have done. Another potential difficulty is 
that the AH location may jump discontinuously, or multiple horizons may form. We believe 
this problem can be handled by simply tracking newly formed horizons and phasing in a new 
AH condition in place of the old one(s). Multiple black holes can be handled by a variation 
of the HLC. We are beginning to examine these issues now with both the 2D axisymmetric 



NCSA black hole code [13|, and a 3D code using harmonic slicing |24], and will report on 



this work in future papers. 
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FIGURES 

FIG. 1. The evolution of the radial metric function A and the coordinate position of the AH are 
plotted for both the BHS code and the HLC code. The dashed line marked "Free Horizon" shows 
the coordinate position of the AH without the HLC, while the dashed line is the "Locked Horizon". 
The solid lines marked "Free A" show the evolution (in time intervals of lOM) of the radial metric 
function A with no shift, while the solid lines, which do not change after 5M, represent the locked 
case. 

FIG. 2. (a) For each of the 4 runs we show a measure of the error in the hamiltonian constraint 



as a function of time. Each point represents \/j2''j=T Gtt{jY, i-e., the sum of the error outside 
the horizon. The solid (dashed) lines are obtained with (without) the AH scheme, (b) The mass 
of the AH is shown as a function of time. 



